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Abstract 

Systems of nearly parallel, slender vortex filaments in which angular momentum is conserved 
are an important simplification of the Navier-Stokes equations where turbulence can be studied 
in statistical equilibrium. We study the canonical Gibbs distribution based on the Klein-Majda- 
Damodaran (KMD) (|7j) model and find a divergence in the mean square vortex position from 
that of the point vortex model of at moderate to high temperature. We subsequently develop 
a free-energy equation based on the non-interacting case, with a spherical constraint, which we ap- 
proximate using the method of Kac-Berlin adding a mean-field term for logarithmic interaction. 
We use this free-energy equation to predict the Monte Carlo results. 
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I. INTRODUCTION 



In statistical equilibrium the behavior of collections of nearly parallel slender vortex 
filaments, periodic in the z-direction and confined by angular momentum in the xy-directions, 
is understood to be nearly identical to that of point vortices (in which no z-direction exists) 
of Onsager for a wide range of temperatures 16|. When a filament is quite straight, 
its internal configuration has a much smaller influence on macroscopic statistical properties 
of the system as a whole (e.g. density and shape of a bundle of filaments) than the mean 
position of its center-line relative to the center-lines of the other filaments in the bundle. 
Taking advantage of this fact, j^] were able to derive a mean- field theory for the density 
of a system of N nearly parallel vortex filaments using the asymptotically derived PDE 
of Klein-Majda-Damodaran (KMD) 3|,|l4], adding mathematical rigor to the statistical 
equilibrium models for turbulence that Chorin has pioneered ^|. However, no adequate 
theory yet exists for the behavior of the KMD system at moderate-to-high temperature 
where turbulent behavior is more extreme and long-range order is minimal. Therefore, we 
must rely on computational methods to explore this important regime. Using the Path 
Integral Monte Carlo method of jj, originally created to study quantum bosons, we sample 
the Gibbs canonical ensemble for the Hamiltonian that []| have derived. Our findings show 
that the point vortex analogy indeed breaks down at high temperatures. While the density 
of point vortex systems increases monotonically with increasing temperature, we show that 
the density of the three-dimensional system actually decreases at a particular temperature. 

As Nordborg and Blatter point out there is a close relationship between the quan- 
tum mechanics of bosons in (2+l)-D and the behavior of vortex filaments. Using methods 
developed for quantum field theory to evaluate functional integrals, we derive a mean-field 
free-energy valid for all temperatures that predicts the behavior of filaments represented by 
smooth curves (a continuum model). 
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II. MODEL 



A. Slender Vortex Filaments 

A system of N Nonlinear Schroedinger Equations (NLSEs) describes the time-evolution 
of slender vortex filaments 

-i^.a^. + l^l^, (1) 

where ipj(a,t) = Xj(a,t) + iyj(a,t) is the position of vortex j at position a along its length 
at time t, and a is the core structure constant (0). Vortex strengths are assumed to all be 
the same and are set to unity. The position in the complex plane, ipj(u, t), is assumed to be 
periodic in a with period L. 

This system of PDEs can be expressed as a Hamiltonian system 

r L N 1 f) I ( \ 2 f L N N 

Hn = oc Y.^-^ d(J - EE^I^W - M<r)\*r- 
Jo k=i 1 ° a Jo k=i i>k 

To this we also add a trapping potential which conserves angular momentum, 

In= / ^|^.(a)| 2 rfa. (2) 
Jo k=i 

For our simulations we assume that the filaments are piecewise linear, divided into an 
equal number of segments of equal length. This discretization leads to the Hamiltonian, 



M N | , . , / / -M2 M N N 

h n( m) =a j:j: i '^ + 1 >-^>i -EEE* lo s 1*0) " 

j=l k=l j=l k=l i>k 

and angular momentum 

M N 
3=1 k=l 

where 5 is the length of each segment and M is the number of segments. For purposes of 
later discussion, the point where two segments meet is called a "bead" in PIMC terminology. 
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B. Path Integrals and Partition Functions 



The path integral method of Feynman in its imaginary time density matrix format in- 
volves evaluating the Gibbs measure of a set of paths, 

Gn[M) ~ Zn(M) ' (4) 

where 

Z N (M)= J2 °n{M). (5) 

allpath^ 

The Gibbs probability, used in the p| model by [9], gives a probability for a path, and that 
allows us to use the path integral method. Originally developed for quantum systems of 
bosons in imaginary time, PIMC applies to the vortex model perfectly as shown by [lol ]. 

Although we use most of Ceperley's original PIMC methods, there is at least one major 
difference in notation. In quantum path integral computations (3 becomes the imaginary 
time length of the path. Since we are modeling real filaments with their own periodic 
length, L, it is important to understand that our (3 corresponds to 1/h and not time. 

III. NUMERICAL RESULTS 

This section describes our numerical results based on the model and method described 
above. Filaments were initialized by scattering them with uniform randomness in a square of 
side 10. Monte Carlo moves involved choosing first a vortex filament to change with uniform 
randomness, then choosing the type of move, either moving the entire chain or rearranging 
the internal configuration via bisection. If moving the entire chain, a new point was chosen 
within a square with a side-length of 10, centered on at the current filament's xy-planar 
position. 

nergy was calculated the same way for both types of moves, using the multilevel method 
Therefore, even the wholechain move had the possibility of being rejected before 
energy was fully evaluated. 

We ran our Monte Carlo simulations until energy settled down to a steady mean. For 
the high (3 (i.e. low temperature) simulations, we see triangular lattices form upon conver- 
gence. Although the fluid remains in a liquid state, these lattices have a crystaline structure 
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resembling that of a solid in which the filaments vibrate but maintain a fixed position w.r.t. 
their neighbors (Figure ^ right). 

Our first finding is a confirmation of findings in ^ which demonstrate the relationship 
between square containment radius, R = (max^l^l 2 )), and the parameters (3 and \i and the 
circulation Q, where a was allowed to remain fixed. We find near perfect agreement in the 
line slopes to formula for square containment radius R 2 = Nf3/(47t{i). Additionally, the 
Hamiltonian used in is divided by ir, which is not done in [jj. Therefore, the formula for 
our square containment radius is 



R 2 = N(3/(4fi), 



(6) 



which gives the slope equal to /3/(4/x) in agreement Figured] left. 
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FIG. 1: Low-temperture ensures a well-defined containment radius. 

Our next finding, shown in Figure El is much more important as it deviates from the 
formula of ^ . Here we see that the pure interaction of the trapping potential term and the 
logarithmic interaction term is no longer valid. While the asymptotic assumptions of the 
model are not broken, 

L > A 2 > e, (7) 

where A (Figure EJ) is the amplitude of the vortices and e is the core size (here assumed to 
be small), the mean-field behavior at high (3 breaks down. At the point where the slope of 
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the curve changes around (3 ~ 0.004, the entropy begins to affect the statistical equilibrium 
R 2 of the bundle, reducing the slope. As we will see in Section 4, the discretization has a 
profound lowering effect on this part of the curve, making a continuum model necessary for 
accurate prediction of the entropy regime f3 < 0.004. 

Mean Square Vortex Position Vs.p 
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FIG. 2: At a specific temperature {(3 ~ 0.004) the entropy takes over from logarithmic interaction 
as the main force for expansion. Here L = 100, a = 10 6 , fj, = 2000, N = 200, and M = 64. 

IV. FREE ENERGY THEORY 

While the above results are interesting of themselves, our investigation is incomplete with- 
out a free-energy theory to explain them for all choices of parameters. Directly calculating 
the free-energy of the interacting system is impossible with current mathematical knowledge. 
However, we might approximate the behavior of R 2 under a mean-field interaction term if 
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FIG. 3: The amplitude is much less than the length L = 100. 

we place the filaments in the non-interacting system on average at an undetermined distance 
R from the center via a spherical constraint. Quantum mechanical path integral methods 
are particularily useful for this section. 

The classical mean-field action (in imaginary time causing the potential to flip sign) is 

S mf = £dr - ^a(3\dMr)\ 2 ~ /i|^(r)| 2 + /3iVlog(it: 2 )/4, (8) 
and the spherical constraint, we enforce with a delta function 

W(r)| 2 -n (9) 

which has integral representation 

/CO J 
^exp-^(T)| 2 -n (10) 
■CO ^ 
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at each r. 
Now let 

gsphere = ^ ^ _ 1 [ a/W ( T )f + + 2/U )|^(t)| 2 - 2/3iV log^ 2 )^ - ^V(t)] , (11) 

be the action. 

The non-dimensional free-energy (/ = j3F, where F is in energy units) as a function of 
ia{r) is then 

/H = -ln J D^ex V S s $ ere , (12) 

where Dip represents the functional integral over all paths ip{r) and the partition function 
for N filaments is 

Z N = J Daexp(-Nf[ia}) (13) 

j^j]. (The exchange of the functional integration of cr with that of ip is permitted in this case 
because the action is negative definite.) Clearly as N — > oo the saddle point gives the main 
contribution [lj. Therefore, 

/«,= lim -l\nZ N = f[r]}, (14) 

N^oo iv 

where r] = ia(r) is the saddle point. We let f3' = f3N and a' = aN^ 1 to take the non- 
extensive limit while keeping a(3 = a'/3' constant. 

The free-energy involves a simple harmonic oscillator with a constant external force 

f[ia) = - In h[ia) + I dr - \%oR 2 - (3' log(i? 2 )/4. (15) 
Jo 2 

Here h is the partition function for a quantum harmonic oscillator in imaginary time, 

h[ia]= y^exp^ dr - ^m[\d T ip\ 2 + cu 2 \iP\ 2 }, (16) 

which has the well known solution for periodic paths in (2+l)-D where we have integrated 
end-points over the whole plane, 

where m = a' f3' and u 2 = {ia + 2/x)/(a'^') Q, Q- 



S 



Let us first make a change of variables A = ia + 2/x. Then the free-energy reads 



/[A] 



rfr(/i - -\)R 2 



- (3'L\og(R 2 )/A-\n- 



-ujL 



-ujL 



-i) 



2 ' 



(18) 



where u; = \/A/ TafW). 

n 

Letting 77 be the saddle point and independent of r |5| , we find that the free-energy is 



M = - \v)LR 2 - P'L log( J R 2 )/4 + JJL + 2 log 

where df/dX = and df/dR 2 = at 77 and i? 2 . 
Since 

5/ 



exp 



a'(5' 



£1-1 



( /J -A/2)L-L/3V(4/? 2 ), 



(19) 



(20) 



I? 2 



4^ -A/2)' 



Substituting for R 2 we get 



/[A]=/5'L/4 + 



A 



Taking the derivative 



L + 2 lo£ 



A 



exp 



:L -1 



log 



L 



df 



4 °4(/i-A/2) 



LP' 



exp 



K/i-A/2) : 



(21) 



(22) 



(23) 



we get a transcendental equation that cannot be solved analytically. However, it is clear 
that A = 2/i and A = are transition points. For the rest of the paper we will drop the 
primes from (3' and a'. 

This system is a quantum one in all but name with an inverse quantum temperature of 
L. To study any phase transitions, we must take the L — > 00 limit because quantum phase 
transitions only occur at absolute zero Q|. Therefore, taking the limit we get a per unit 
length energy of 

/r-W = \ + V^im - f log (45^)) • (24) 

where 

77 = 2/x - i/3(-/3 2 a ± v 7 /? 4 " 2 + 32a/3/j), (25) 
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of which we take 

77 = 2/i— \(3{-(3 2 a + + (26) 

8 

as giving physical results (shown in Figure 0) 

The two points where the free-energy is non-analytic, 

V = 0, 2fi, (27) 

correspond to 

e = -ffi 0. <») 

which we will call /5_ c and /3 C respectively 

The change in the behavior of R 2 we observe in our Monte Carlo results in Figure El at 
(3 ~ 10~ 2 is not a phase transition but a smooth transition, since the free-energy is smooth 
at all (3 other than /3_ c and j3 c . 

The specific heat 

C = -/? 2 ^? (29) 
[l| (see Figure HJ) can be represented in series as 



C = --v^/i/a/T 1 / 2 + O(higher), (30) 

indicating a critical exponent of v = —1/2, suggesting a second-order (continuous) phase 
transition. Ill 



V. CONCLUSION 

To our knowledge no one has done Monte Carlo simulations of the system proposed in (t| . 
Excellent simulations have been done in cases where boundaries are periodic in all directions 



using the London energy functional for flux line lattices, which differs 

n 



such as 10] and Ibj 

n 

from that of [7j only in that the interaction potential is a modifed Bessel's function (log-like at 
short distances). However, free boundary conditions with the addition of the conservation of 
angular momentum make this problem different and specifically applicable to fluid statistics 
and to a continuous filament model of the Gross-Petaevskii equation in the Thomas limit - 
details of which will appear in a later paper. 

Our findings - comparing the point vortex expression for R 2 with|2I]- indicate a special 
regime of temperature where entropy of the filaments has a profound effect on variance 
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FIG. 4: The specific heat in Equation 1291 shown here with parameters a = 10 6 , N = 200, and 
fj, = 2000 diverges at /3 = and (3 = — \j For (3 in this interval it has complex value. 



of vortex position, which has not been observed before. Also new is the discovery of a 
quantum, i.e. infinite L, phase transition in the continuous vortex filaments model at zero 
value of inverse temperature as well as a negative (3 transition in Equation EHJ According to 
Sachdev characteristics of the quantum transition extend to finite values of L. Hence, 
our comparison between theoretical values for R 2 at infinite L and PIMC values of R 2 at 
finite L are valid. 

We developed a theory based on a combination of the non-interacting case with a mean- 
field estimate for logarithmic interaction. The resulting free-energy allows us to predict the 
(3 point where entropy forces R 2 to increase with temperature and becomes the dominant 
resistance to the angular momentum rather than interaction. We have also shown that this 
is not a phase transition for this model. More interesting behavior may be observed at 
negative (3. We also plan to apply our methods to periodic boundaries. 



Lions and Majda 



fl 



have already suggested applicability of their novel derivations in the 



area of geophysical and astrophysical convert ion such as p] have modeled. We say thai our 
results are equally applicable and we may extend them to Bose-Einstein Condensates in the 
future. 
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FIG. 5: As the broken segment model used in the Monte Carlo simulations approaches the contiuum 
limit M = oo, the change in R 2 with f3 is predicted to approach a "v" -shaped (corresponding to a 
polynomial function in 0) curve. 
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